Nesta lista, vamos trabalhar com o seguinte processo estocástico AR(1):
\[ z_t = \rho z_{t-1} + \varepsilon_t , \]
com \(\varepsilon_t \sim N(\mu,\sigma^2)\). Inicialmente, assumiremos que \(\mu=0\), \(\rho = 0.95\) e \(\sigma = 0.007\), calibração de Cooley & Prescott (1995).
Ao discretizar um processo que assume valores contínuos, é necessário determinar o intervalo do grid, o número de pontos, e sua localização. Além disso, é necessário determinar probabiliades condicionais de transição entre os estados de forma consistente com o processo original.
Pelo método de Tauchen, o intervalo de um grid de \(N\) pontos é determinado por \([\mu - \theta_N, \mu + \theta_N]\), em que \(\theta_N = m \frac{\sigma}{\sqrt{1-\rho^2}}\) é o desvio padrão do processo estocástico, multiplicado por um parâmetro de escala \(m\). Os demais \(N-2\) pontos são distribuídos uniformemente no intervalo.
N = 9
rho = 0.95
sigma = 0.007
mu = 0
m = 1
thetaN = m * sigma / sqrt(1-rho^2)
thetaT = seq(-thetaN, thetaN, length.out = N)
as.matrix(thetaT)## [,1]
## [1,] -0.022417942
## [2,] -0.016813456
## [3,] -0.011208971
## [4,] -0.005604485
## [5,] 0.000000000
## [6,] 0.005604485
## [7,] 0.011208971
## [8,] 0.016813456
## [9,] 0.022417942Seja \(\theta_{i,t}\) um estado no período \(t\), e \(\theta_{j, t+1}\) um estado no período \(t+1\). No método de Tauchen, as probabilidades de transição dos pontos interiores do grid são dadas por
\[ p_{i,j} = \Phi \left( \frac{\theta_j + \Delta \theta/2 - (1-\rho) \mu - \rho \theta_i}{\sigma} \right) - \Phi \left(\frac{\theta_j - \Delta \theta/2 - (1-\rho) \mu - \rho \theta_i}{\sigma}\right)\]
em que \(\Phi\) é a CDF da \(N(\mu, \sigma^2)\). Nos cantos do grid, temos que
\[ p_{i,1} = \Phi \left(\frac{\theta_1 - (1-\rho) \mu - \rho \theta_i + \Delta \theta/2}{\sigma} \right) \]
\[ p_{i,N} = 1 - \Phi \left(\frac{\theta_N - (1-\rho) \mu - \rho \theta_i - \Delta \theta/2}{\sigma} \right).\]
delta = 2 * thetaN / (N-1)
PT = matrix(NA_real_, 9, 9)
for(i in 1:N) {
for(j in 1:N) {
if(j == 1) {
PT[i,j] = pnorm( (thetaT[1] - (1-rho)*mu - rho*thetaT[i] + delta/2) / sigma )
} else if(j==N) {
PT[i,j] = 1 - pnorm( (thetaT[N] - (1-rho)*mu - rho*thetaT[i] - delta/2) / sigma )
} else {
PT[i,j] =
pnorm((thetaT[j] + delta/2 - (1-rho)*mu - rho*thetaT[i]) / sigma) -
pnorm((thetaT[j] - delta/2 - (1-rho)*mu - rho*thetaT[i]) / sigma)
}
}
}
round(PT, 4)## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9]
## [1,] 0.5949 0.2561 0.1162 0.0287 0.0038 0.0003 0.0000 0.0000 0.0000
## [2,] 0.3014 0.3090 0.2498 0.1099 0.0263 0.0034 0.0002 0.0000 0.0000
## [3,] 0.1001 0.2154 0.3101 0.2432 0.1038 0.0241 0.0030 0.0002 0.0000
## [4,] 0.0206 0.0867 0.2225 0.3108 0.2365 0.0979 0.0220 0.0027 0.0002
## [5,] 0.0025 0.0201 0.0922 0.2296 0.3111 0.2296 0.0922 0.0201 0.0025
## [6,] 0.0002 0.0027 0.0220 0.0979 0.2365 0.3108 0.2225 0.0867 0.0206
## [7,] 0.0000 0.0002 0.0030 0.0241 0.1038 0.2432 0.3101 0.2154 0.1001
## [8,] 0.0000 0.0000 0.0002 0.0034 0.0263 0.1099 0.2498 0.3090 0.3014
## [9,] 0.0000 0.0000 0.0000 0.0003 0.0038 0.0287 0.1162 0.2561 0.5949rowSums(PT) # Linhas da matriz devem somar 1.## [1] 1 1 1 1 1 1 1 1 1A definição do grid no método de Rouwenhorst é um caso particular do método de Tauchen, com \(m = \sqrt{N-1}\).
thetaN = sqrt(N-1) * sigma / sqrt(1-rho^2)
thetaR = seq(-thetaN, thetaN, length.out = N)
as.matrix(thetaR)## [,1]
## [1,] -0.06340751
## [2,] -0.04755564
## [3,] -0.03170376
## [4,] -0.01585188
## [5,] 0.00000000
## [6,] 0.01585188
## [7,] 0.03170376
## [8,] 0.04755564
## [9,] 0.06340751A matriz de transição \(P\), porém, é obtida recursivamente.
p = (1+rho)/2
P_init = matrix(c(p, 1-p, 1-p, p), 2, 2, byrow = TRUE)
for(i in 3:N) {
PR =
rbind(cbind(p*P_init, 0), 0) +
rbind(cbind(0, (1-p)*P_init), 0) +
rbind(0, cbind((1-p)*P_init, 0)) +
rbind(0, cbind(0, p*P_init))
PR = PR / rowSums(PR)
P_init = PR
}
round(PR, 4)## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9]
## [1,] 0.8167 0.1675 0.0150 0.0008 0.0000 0.0000 0.0000 0.0000 0.0000
## [2,] 0.0209 0.8204 0.1469 0.0113 0.0005 0.0000 0.0000 0.0000 0.0000
## [3,] 0.0005 0.0420 0.8231 0.1261 0.0081 0.0003 0.0000 0.0000 0.0000
## [4,] 0.0000 0.0016 0.0630 0.8247 0.1051 0.0054 0.0001 0.0000 0.0000
## [5,] 0.0000 0.0001 0.0032 0.0841 0.8253 0.0841 0.0032 0.0001 0.0000
## [6,] 0.0000 0.0000 0.0001 0.0054 0.1051 0.8247 0.0630 0.0016 0.0000
## [7,] 0.0000 0.0000 0.0000 0.0003 0.0081 0.1261 0.8231 0.0420 0.0005
## [8,] 0.0000 0.0000 0.0000 0.0000 0.0005 0.0113 0.1469 0.8204 0.0209
## [9,] 0.0000 0.0000 0.0000 0.0000 0.0000 0.0008 0.0150 0.1675 0.8167rowSums(PR) # Linhas da matriz devem somar 1.## [1] 1 1 1 1 1 1 1 1 1set.seed(96452)
t_span = 10000
eps = rnorm(t_span, mu, sigma)
z = rep(0, t_span)
# Processo contínuo
zc = z
for(t in 2:t_span)
zc[t] = rho * zc[t-1] + eps[t]
# Método de Tauchen
zt = z
for(t in 2:t_span) {
i = which(thetaT == zt[t-1])
j = sum(cumsum(PT[i, ]) < pnorm(eps[t], mu, sigma)) + 1
zt[t] = thetaT[j]
}
# Método de Rouwenhorst
zr = z
for(t in 2:t_span) {
i = which(thetaR == zr[t-1])
j = sum(cumsum(PR[i, ]) < pnorm(eps[t], mu, sigma)) + 1
zr[t] = thetaR[j]
}Os gráficos a seguir exibem o processo contínuo e as discretizações de Tauchen e de Rouwenhorst para os cem primeiros períodos:
Aparentemente, a discretização de Tauchen funciona melhor com os parâmetros dados. Os pontos mais próximos no grid de Tauchen capturam de forma mais fina o movimento do processo, sem perder muitos pontos extremos. O método de Rouwenhorst funciona melhor para valores extremos.
Com 25 pontos, o método de Rouwenhorst produz resultados semelhantes ao de Tauchen.
fit_zt = lm(zt ~ shift(zt) + 0)
fit_zr = lm(zr ~ shift(zr) + 0)
fit_zrl = lm(zrl ~ shift(zrl) + 0)
coefs = c(fit_zt$coefficients,
fit_zr$coefficients,
fit_zrl$coefficients)
sigmas = c(summary(fit_zt)$sigma,
summary(fit_zr)$sigma,
summary(fit_zrl)$sigma)
tbl = data.frame(method = c("Tauchen 9 pts",
"Rouwenhorst 9 pts",
"Rouwenhorst 25 pts"),
rho = coefs, sigma = sigmas)
tblPelo método de Tauchen com 9 pontos no grid, estimamos \(\hat{\rho} = 0.8891\), valor distante do \(\rho\) verdadeiro. Aumentando o número de pontos no grid para
thetaN = m * sigma / sqrt(1-rho^2)
thetaT = seq(-thetaN, thetaN, length.out = 25)
delta = 2 * thetaN / (25-1)
PT = matrix(NA_real_, 25, 25)
for(i in 1:25) {
for(j in 1:25) {
if(j == 1) {
PT[i,j] = pnorm( (thetaT[1] - (1-rho)*mu - rho*thetaT[i] + delta/2) / sigma )
} else if(j==25) {
PT[i,j] = 1 - pnorm( (thetaT[25] - (1-rho)*mu - rho*thetaT[i] - delta/2) / sigma )
} else {
PT[i,j] =
pnorm((thetaT[j] + delta/2 - (1-rho)*mu - rho*thetaT[i]) / sigma) -
pnorm((thetaT[j] - delta/2 - (1-rho)*mu - rho*thetaT[i]) / sigma)
}
}
}
# Método de Tauchen
ztl = z
for(t in 2:t_span) {
i = which(thetaT == ztl[t-1])
j = sum(cumsum(PT[i, ]) < pnorm(eps[t], mu, sigma)) + 1
ztl[t] = thetaT[j]
}
fit_ztl = lm(ztl ~ shift(ztl) + 0)
rbind(tbl, cbind(method = "Tauchen 25 pts", rho = fit_ztl$coefficients, sigma = summary(fit_ztl)$sigma))